arg = argv;
X = load("-ascii", arg{1})';
x = X(1,:);
sample = X(7,:);

% smooth, subtract baseline, and calculate the derivatives
order = 2;
signal = dwfilt(sample, order); # Savitsky-Golay with Durbin-Watson criterion, order = 2

% find minima, skipping every other one, then filter the resulting curve
#mix = findminima(signal);
#baseline = interp1(mix, signal(mix), 1:length(signal));
mix = findminima(signal);
baseline = 0.5*sgolayfilt(interp1(mix, signal(mix), 1:length(signal)), 2, 19);
#signal = signal - baseline;
d1 = dwfilt(diff(signal), 1);
d2 = dwfilt(diff(d1), order);
d3 = dwfilt(diff(d2), order);

% trim all signals to the length of the 3rd derivative
x = x(3:length(d3)+2);
baseline = baseline(3:length(d3)+2);
signal = signal(3:length(d3)+2);
d1 = d1(2:length(d3)+1);
d2 = d2(2:length(d3)+1);

% set thresholds
snoise = noise(sample);
thr00 = 2*snoise; # the lower limit for height
thr0 = 5*snoise;
thr1 = 3*noise(d1);
thr2 = 3*noise(d2);

% find the negative troughs of the 2nd derivative
region = find(d2 < -thr2 & signal > thr0 & (d1 > thr1 | d1 < -thr1) );
mark = zeros(1, length(x));
mark(region) = ones(1, length(region));

% bridge the gaps (if two marked regions are one point apart
% it's likely because d1 constraint above knocked out a point near zero)
for i = 2:length(mark)-1
  if (mark(i) == 0 & mark(i-1) == 1 & mark(i+1) == 1)
    mark(i) = 1;
  endif
endfor

starts = find(diff(mark) > 0);
ends = find(diff(mark) < 0);
% Kludge !!!
if ends(1) < starts(1)
  if ends(2) < starts(2)
    ends = ends(2:end);
  else
    disp("bad segment order");
    exit(-1);
  endif
endif
common_size = min(length(starts), length(ends));
starts = starts(1:common_size);
ends = ends(1:common_size);

% extend the peak regions by annexing the areas of substantial 1st derivative
region = {}; # reuse region array in a different way
for i = 1:length(starts)
  from = starts(i);
  to = ends(i);
  if to - from <= 3 # discard short regions as noise
    continue;
  endif

  fusedbefore = fusedafter = 0;

  j = from;
  while ( d1(j) > thr1 & signal(j) > thr0 )
    if (j == 1)
      break;
    endif
    if ( !fusedbefore & sign(d3(j)) != sign(d3(j-1)) )
      fusedbefore = j;
    endif
    j--;
  endwhile
  from = j;
  % don't bother with small fused areas
  if fusedbefore - from <= 7
    fusedbefore = 0;
  endif

  j = to;
  while ( d1(j) < -thr1 & signal(j) > thr0 )
    if ( !fusedafter & sign(d3(j)) != sign(d3(j-1)) )
      fusedafter = j;
    endif
    j++;
    if j > length(signal) 
      break;
    endif
  endwhile
  to = j-1;
  if to - fusedafter <= 7
    fusedafter = 0;
  endif

  if ( fusedbefore )
    region{end+1} = [from, fusedbefore-1];
    from = fusedbefore;
  endif
  if ( fusedafter )
    region{end+1} = [from, fusedafter];
    region{end+1} = [fusedafter+1, to];
  else
    if to - from > 5
      region{end+1} = [from, to];
   endif
 endif
endfor

% this detection routine somehow resuls in overlapping
% regions. Join those that overlap by more than one point
if length(region) > 1
  region1 = region;
  region = {};
  i = 1;
  while ( i < length(region1) )
    reg1 = region1{i}(1):region1{i}(2);
    reg2 = region1{i+1}(1):region1{i+1}(2);
    newreg = intersect(reg1, reg2);
    if length(newreg) >= 2
      newreg = union(reg1, reg2);
      region{end+1} = [newreg(1), newreg(end)];
      i = i + 2;
    else
      region{end+1} = [reg1(1), reg1(end)];
      i++;
    endif
  endwhile
  if length(intersect(region{end}(1):region{end}(2), region1{end}(1):region1{end}(2))) <= 1
    region{end+1} = region1{end};
  endif
endif

% calculate and save peak data
count = 0;
file = fopen(sprintf("%s.peak", arg{2}), "w");
abifile = fopen(sprintf("%s.abi", arg{2}), "w");
format short g;
for i = 1:length(region)
  range = region{i};
  from = range(1);
  to = range(2);
  before = max(1, from - 2);
  after = min(length(x), to + 2);
  [mind2, mind2x] = min(d2(from:to));
  sizex = from + mind2x - 1;
  if mind2 > -thr2 # this is not a real peak
    sizex = floor((from + to)/2);
  else
    % stop toward real maximum if there is an offest in the 2nd derivative
    if sizex != from & sizex != to
      if signal(sizex - 1) > signal(sizex)
	sizex = sizex - 1;
      endif
      if signal(sizex + 1) > signal(sizex)
	sizex = sizex + 1;
      endif
    endif
  endif

% this code calculates artificial baseline
  edge1 = signal(from)/2;
  edge2 = signal(to)/2;
  base = edge1 + (sizex - from)*(edge2 - edge1)/(to - from);
  height = signal(sizex) - base;
  if height < thr00 # height small or negative? go away.
    continue;
  endif
  basearea = mean([edge1, edge2])*(to - from);
  area = sum(signal(from:to)) - basearea;

#   % this works with the baseline estimated from data
#   edge1 = signal(from) - baseline(from);
#   edge2 = signal(to) - baseline(from);
#   base = edge1 + (sizex - from)*(edge2 - edge1)/(to - from);
#   height = signal(sizex) - base;
#   if height < thr00 # height small or negative? go away.
#     continue;
#   endif
#   basearea = mean([edge1, edge2])*(to - from);
#   area = sum(signal(from:to) - baseline(from:to)) - basearea;

  if x(sizex) < 30 
    continue;
  endif
  count++;
  fprintf(file, "%d\t%d\t%d\t%6.2f\t%6.2f\t%6.2f\n", from, to, sizex, x(sizex), height, area);
  fprintf(abifile, "B,%d\t%d\t%6.2f\t%6.2f\t%6.2f\t%d\n", count, sizex, x(sizex), height, area, sizex);
endfor
fclose(file);
fclose(abifile);
printf("%d peaks", count);
  
X = vertcat(x, 1 + 0*baseline, signal, d1, d2, d3)';
save("-ascii", sprintf("%s.deriv", arg{2}), "X");

